Code
source("requirements.R")Atelier Méthodes CED
Ce document détaille pas à pas les étapes de réalisation des analyses géométriques les plus simples quand on débute dans l’application de cette méthode avec R studio :
Les exemples que nous utiliserons ici sont tirés de la littérature académique (avec accès aux données utilisées par leurs auteur/ices) :
NB : des éléments de syntaxe de R peuvent être retrouvés ici : https://iqss.github.io/dss-workshops/R/Rintro/base-r-cheat-sheet.pdf
Dans notre cas, il est possible de simplifier l’installation et le chargement des packages en faisant appel au fichier requirements.R qui liste l’ensemble des éléments nécessaires. (À condition qu’il soit bien dans le dossier depuis lequel vous exécutez votre code R)
source("requirements.R")Exécuter la commande ci-dessus revient à installer les packages et les charger, processus sinon détaillé ci-dessous.
script pour chaque package : install.packages(“NomPackage”)
Liste des packages les plus utiles pour ce qui nous concerne :
# install.packages("Factoshiny")
# install.packages("FactoMineR")
# install.packages("shiny")
# install.packages("ggplot2")
# install.packages("FactoInvestigate")
# install.packages("ggpp")
# install.packages("dplyr")
# install.packages("ggpmisc")
# install.packages("ggrepel")
# install.packages("readxl")
# install.packages("haven")
# install.packages("GDAtools")
# install.packages("sf")
# install.packages("questionr")script pour chaque package : library(“NomPackage”)
# library(Factoshiny)
# library(FactoMineR)
# library(shiny)
# library(ggplot2)
# library(FactoInvestigate)
# library(ggpp)
# library(dplyr)
# library(ggpmisc)
# library(ggrepel)
# library(readxl)
# library(haven)
# library(GDAtools)
# library(sf)
# library(questionr)
# TODO: vérif mais installé et chargé par défaut non ?
# Peut virer pour nettoyer ?
#library(stats)
#library(base)Un mini jeu de données à 2 variables continues
–> Taux de chômage en France 2000 / 2007 (données OCDE)
Opération possible via le menu déroulant de Rstudio (le fichier est dans cet exemple au format excel) :
File –> Import Dataset –> From Excel
Ou sinon directement en code en exécutant la commande suivante :
# ici il faut possiblement adapter le chemin à sa propre arborescence
ex1 <- read_excel("data/EX1_DonneesOCDE_Base2.xls")ex1 = c’est le nom choisi ici qui servira à appeler le jeu de données dans R dans la suite de cet exemple. Vous pouvez tout aussi bien l’appeler ‘fraise’ si vous le souhaitez – vous l’écrivez alors ainsi dans le script :
fraise <- read_excel("ex5b-trivial-acp.xls")
NB : Il est plus simple de choisir un nom court, facile à manipuler.
Je peux ensuite visualiser mon fichier avec View() ou afficher les premiers éléments avec head() :
# View(ex1) # Décommenter cette ligne pour utiliser view
# view vous sera utile dans Rstudio ou équivalent
# mais comme il ouvre un autre objet, on se contente ici d'afficher
# les premières lignes du fichier avec head()
head(ex1)# A tibble: 6 × 5
id TxChmge2000 TxChmge2007 TxChmge2000CR TxChmge2007CR
<chr> <dbl> <dbl> <dbl> <dbl>
1 Paris 10.2 9.16 0.126 0.871
2 Seine-et-Marne 6.55 6.64 -1.09 -0.556
3 Yvelines 6.31 6.38 -1.16 -0.707
4 Essonne 6.45 6.12 -1.12 -0.855
5 Hauts-de-Seine 7.93 7.55 -0.632 -0.0422
6 Seine-Saint-Denis 13.2 11.5 1.11 2.20
Je peux également afficher les principales stats (min, max, médiane, moyenne, quartiles) avec summary()
summary(ex1) id TxChmge2000 TxChmge2007 TxChmge2000CR
Length:95 Min. : 5.034 Min. : 4.478 Min. :-1.5852
Class :character 1st Qu.: 8.069 1st Qu.: 6.400 1st Qu.:-0.5856
Mode :character Median : 9.277 Median : 7.341 Median :-0.1879
Mean : 9.848 Mean : 7.625 Mean : 0.0000
3rd Qu.:10.825 3rd Qu.: 8.556 3rd Qu.: 0.3220
Max. :21.472 Max. :12.205 Max. : 3.8279
TxChmge2007CR
Min. :-1.7831
1st Qu.:-0.6942
Median :-0.1611
Mean : 0.0000
3rd Qu.: 0.5272
Max. : 2.5944
ggplot(ex1, aes(x = TxChmge2000CR, y = TxChmge2007CR)) +
geom_point(color = "#EF6A40", size = 3) +
geom_text(aes(label = id), hjust = -0.1, vjust = 0.5, size = 3) +
theme_minimal() +
labs(
title = "Ici je peux mettre un titre",
x = "légende pour l'axe x",
y = "légende pour l'axe y"
)Dico script R
Trucs et astuces
ggplot(ex1, aes(x = TxChmge2000CR, y = TxChmge2007CR)) +
geom_point(color = "steelblue", size = 3) +
geom_text_repel(aes(label = id), size = 3, max.overlaps = Inf) +
theme_minimal() +
labs(
title = "Ici je peux mettre un titre",
x = "légende pour l'axe x",
y = "légende pour l'axe y"
)==> J’appelle ce graphe intermédiaire “p”
Dans le script ci-dessous, c’est ce que signifie : ==> p <-
# Identifier les points extrêmes à annoter (par exemple sur Y)
points_a_annoter <- ex1 %>%
arrange(TxChmge2007CR) %>%
slice_head(n = 5) %>% # les 5 plus bas
bind_rows(
ex1 %>% arrange(desc(TxChmge2007CR)) %>% slice_head(n = 7) # les 7 plus hauts
)
# Ajouter aussi le point avec la valeur maximale sur l'axe X
points_max_pop <- ex1 %>% arrange(desc(TxChmge2000CR)) %>% slice_head(n = 9)
# Combiner avec les autres extrêmes
points_a_annoter <- ex1 %>%
arrange(TxChmge2000CR) %>%
slice_head(n = 7) %>%
bind_rows(
ex1 %>% arrange(desc(TxChmge2007CR)) %>% slice_head(n = (7)),
points_max_pop
) %>%
distinct() # pour éviter les doublons
# Nuage de points avec étiquettes seulement sur les points extrêmes
p <- ggplot(ex1, aes(x = TxChmge2000CR, y = TxChmge2007CR)) +
geom_point(color = "steelblue", size = 3) +
geom_text_repel(
data = points_a_annoter,
aes(label = id),
size = 3
) +
theme_minimal() +
labs(
title = "Ici je peux mettre un titre",
x = "légende pour l'axe x",
y = "légende pour l'axe y"
)
print(p)library(ggpmisc)
# Nuage de points avec étiquettes seulement sur les points extrêmes
# + ajout de la droite de régression avec geom_smooth
p +
geom_smooth(method = "lm", se = TRUE, color = "darkred") +
stat_poly_eq(
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
formula = y ~ x,
parse = TRUE,
size = 4,
color = "darkred"
) Factoshiny permet de faire une exploration rapide du résultat de l’analyse géométrique.
# script en commentaire pour qu'il ne s'exécute pas dans le notebook
# acp <- Factoshiny(ex1)Qui fournit le script associé aux opérations menées en click-bouton.
On peut ensuite le copier-coller dans la console R : - pour garder la trace de ce que l’on a fait - pour améliorer, via le script, les graphiques initiaux réalisés par Factoshiny
ATTENTION - Pour retourner dans la console R quand on a joué avec Factoshiny (sans que R ne bugge), il ne faut pas oublier de cliquer sur :
On peut ensuite reproduire le code créé par Factoshiny dans R :
res.PCA<-PCA(ex1,quali.sup=c(1),quanti.sup=c(4,5),graph=FALSE)
plot.PCA(res.PCA,choix='var')plot.PCA(res.PCA,invisible=c('ind','ind.sup'),select='contrib 95',cex=0.5,cex.main=0.5,cex.axis=0.5,label =c('quali'))# Affichage des infos résumées
# Ici on éviter pour alléger le notebook
# summary(res.PCA) # Décommenter si besoinComme vu précédemment, il est possible de charger les données via le menu déroulant de Rstudio. Ci-dessous, on passe plutôt par des lignes de code.
safi <- read_excel("data/EX2_ACP_MirnaSAFI_Extract.xls")
head(safi)# A tibble: 6 × 6
app cult nor mix eco inat3
<dbl> <dbl> <dbl> <dbl> <dbl> <chr>
1 0 0.300 0.5 0.667 0.611 Portugal
2 0.5 0.5 0.833 0.833 0.667 Portugal
3 0.5 0.875 1 0.583 0.583 Portugal
4 0.75 0.833 0.333 0.833 0.361 Autres
5 1 0.5 0.583 0.667 0.833 Asie du Sud-Est
6 0 0.5 0.583 0.375 0.417 Portugal
safi = le nom qui servira à appeler le jeu de données dans R dans la suite de l’exercice.
Statistiques descriptives (variables continues) :
# Afficher les infos résumées
# Ici on évite pour alléger le notebook
# summary(safi) # Décommenter si besoinComme vu précédemment, Factoshiny permet de faire une exploration rapide du résultat de l’analyse géométrique.
# script en commentaire pour qu'il ne s'exécute pas dans le notebook
# acp <- Factoshiny(safi)ATTENTION - Penser à cliquer sur “[x] Quit the app” pour retourner dans la console R quand on a joué avec Factoshiny (sans que R ne bugge)
# graphiques initiaux réalisés dans Factoshiny
res.PCA<-PCA(safi,quali.sup=c(6),graph=FALSE)
plot.PCA(res.PCA,choix='var')# Colorier les individus selon le pays
plot.PCA(res.PCA,invisible=c('ind.sup'),habillage=6,cex=0.8,cex.main=0.8,cex.axis=0.8,label =c('quali', cex = 3, col = "black"))# Afficher les barycentres des modalités de variable supplémentaire pays
# (possible de peauffiner pour avoir un résultat plus satisfaisant)
plot.PCA(res.PCA,invisible=c('ind.sup'),cex=1.2,cex.main=1.2,cex.axis=1.2,label =c('quali'))# Afficher les résumés et détails
# Ici on ne lance pas pour alléger le notebook
# summary(res.PCA) # Décommenter si besoin
# dimdesc(res.PCA) # Décommenter si besoinA VOUS DE JOUER !!
On vous remet toutefois ci-dessous une proposition d’exploration telle que vue en séance.
#library(haven)
# bug ici chez viviane si on rappelle pas haven
# Et visiblement chez moi dans Rstudio
# Mais pas dans vscode ???
# TODO: check pourquoi the what of the fuck
# ici haven est utilisé pour lire le fichier .dta
# (déjà chargé plus haut)
epices <- read_dta("data/EX3_Epices.dta")
head(epices)# A tibble: 6 × 26
numquest EPICES EPICESX QUINTEPI Sport_Non Spectacle_Non Vacances_Non
<dbl> <dbl> <dbl> <dbl+lbl> <dbl> <dbl> <dbl>
1 1 21.9 15.9 1 [2.0] 2 1 1
2 2 48.5 44.6 4 [5.0] 1 1 1
3 3 45.6 41.4 4 [5.0] 2 1 2
4 4 48.5 44.6 4 [5.0] 1 1 1
5 5 78.7 77.1 4 [5.0] 2 2 2
6 6 28.4 22.9 2 [3.0] 2 1 1
# ℹ 19 more variables: Contact_Famille_Non <dbl>, Hébergement_Non <dbl>,
# SRevenuUC_Faible <dbl>, Propriétaire_Non <dbl>, SBienRapport_Non <dbl>,
# STerre_Non <dbl>, SActions_Non <dbl>, SEconomies_Non <dbl>,
# STéléphoneMobile_Non <dbl>, SOrdiPortable_Non <dbl>, SVoiture_Non <dbl>,
# SLaveVaisselle_Non <dbl>, SInternet_Non <dbl>, SPatrimoine_Faible <dbl>,
# SEmprunts_Importants <dbl>, Mutuelle_Non <dbl>, Couple_Non <dbl>,
# Difficultés_Financières_Oui <dbl>, Ass_Sociale_Oui <dbl>
# Décommenter la ligne du dessou si l'on veut manipuler dans Factochiny :
# acp <- Factoshiny(epices)# Réaliser l'ACP
res.PCA<-PCA(epices,quanti.sup=c(1,2,3,4,10,12,13,14,15,16,17,18,19,20,21,22),graph=FALSE)
# Sortir une représentation graphique
plot.PCA(res.PCA,choix='var',cex=0.75,cex.main=0.75,cex.axis=0.75,
col.quanti.sup='#0000FF')Warning: ggrepel: 5 unlabeled data points (too many overlaps). Consider
increasing max.overlaps
plot.PCA(res.PCA)# Sortir les stats (ici on éviter pour alléger le notebook) :
# summary(res.PCA) # Décommenter si besoin
# dimdesc(res.PCA) # Décommenter si besoinA VOUS DE JOUER !!
Pour l’ACM, on suggère de faire les graphes dans R-Studio via le package GDAtool.
Quelques ressources :
On propose également ci-dessous une version rapide (et perfectible) du code vu durant la séance.
# Charger les données
duval <- read_excel("data/EX4_Cinema_GDAtools.xlsx")
# Faire l'ACM
mca_duval<-MCA(duval,quanti.sup=c(1),quali.sup=c(17),graph=FALSE)
# Sortir les stats
# summary(mca_duval) # Décommenter si besoin
# dimdesc(mca_duval) # Décommenter si besoin
# Visualiser
plot.MCA(mca_duval, choix='var',col.quali.sup='#006400')plot.MCA(mca_duval,invisible= 'ind',col.quali.sup='#006400',label =c('var','quali.sup'))On peut également aller regarder les eigenvalues et les taux modifiés de Benzecri :
mca_duval$eig eigenvalue percentage of variance cumulative percentage of variance
dim 1 0.21332254 12.307070 12.30707
dim 2 0.18894143 10.900467 23.20754
dim 3 0.10058452 5.802953 29.01049
dim 4 0.09653844 5.569526 34.58002
dim 5 0.08489679 4.897892 39.47791
dim 6 0.08079175 4.661063 44.13897
dim 7 0.07834366 4.519827 48.65880
dim 8 0.07566525 4.365303 53.02410
dim 9 0.06951334 4.010385 57.03448
dim 10 0.06820379 3.934834 60.96932
dim 11 0.06742118 3.889683 64.85900
dim 12 0.06071324 3.502687 68.36169
dim 13 0.05856164 3.378556 71.74024
dim 14 0.05420357 3.127129 74.86737
dim 15 0.05230522 3.017609 77.88498
dim 16 0.04622018 2.666549 80.55153
dim 17 0.04526385 2.611376 83.16291
dim 18 0.04380656 2.527301 85.69021
dim 19 0.04004370 2.310213 88.00042
dim 20 0.03761535 2.170116 90.17054
dim 21 0.03568989 2.059032 92.22957
dim 22 0.03276323 1.890186 94.11976
dim 23 0.02846471 1.642195 95.76195
dim 24 0.02620882 1.512047 97.27400
dim 25 0.02579505 1.488176 98.76218
dim 26 0.02145561 1.237824 100.00000
modif.rate(mca_duval)$raw
eigen rate cum.rate
dim 1 0.21332254 12.307070 12.30707
dim 2 0.18894143 10.900467 23.20754
dim 3 0.10058452 5.802953 29.01049
dim 4 0.09653844 5.569526 34.58002
dim 5 0.08489679 4.897892 39.47791
dim 6 0.08079175 4.661063 44.13897
dim 7 0.07834366 4.519827 48.65880
dim 8 0.07566525 4.365303 53.02410
dim 9 0.06951334 4.010385 57.03448
dim 10 0.06820379 3.934834 60.96932
dim 11 0.06742118 3.889683 64.85900
dim 12 0.06071324 3.502687 68.36169
dim 13 0.05856164 3.378556 71.74024
dim 14 0.05420357 3.127129 74.86737
dim 15 0.05230522 3.017609 77.88498
dim 16 0.04622018 2.666549 80.55153
dim 17 0.04526385 2.611376 83.16291
dim 18 0.04380656 2.527301 85.69021
dim 19 0.04004370 2.310213 88.00042
dim 20 0.03761535 2.170116 90.17054
dim 21 0.03568989 2.059032 92.22957
dim 22 0.03276323 1.890186 94.11976
dim 23 0.02846471 1.642195 95.76195
dim 24 0.02620882 1.512047 97.27400
dim 25 0.02579505 1.488176 98.76218
dim 26 0.02145561 1.237824 100.00000
$modif
mrate cum.mrate
dim 1 54.780529523 54.78053
dim 2 38.080354715 92.86088
dim 3 2.930110291 95.79099
dim 4 2.272738011 98.06373
dim 5 0.846460033 98.91019
dim 6 0.508170692 99.41836
dim 7 0.347287703 99.76565
dim 8 0.206241502 99.97189
dim 9 0.020639644 99.99253
dim 10 0.006017926 99.99855
dim 11 0.001449962 100.00000
Ou encore aller vérifier les descriptions et contributions aux dimensions
# # Décommenter si vous voulez la sortie
# # ici on évite pour alléger le notebook
# dimdescr(mca_duval, vars = NULL, dim = c(1,2),
# limit = NULL, correlation = "pearson",
# na.rm.cat = FALSE, na.value.cat = "NA", na.rm.cont = FALSE,
# nperm = NULL, distrib = "asympt",
# shortlabs = TRUE)
# dimcontrib(mca_duval, dim = c(1,2), best = TRUE)# ici dans le plan 1-2 ==> axes = c(1,2)
ggcloud_indiv(mca_duval, type = "i", points = "all", axes = c(1,2),
col = "dodgerblue4", point.size = 0.5, alpha = 0.6,
repel = FALSE, text.size = 2,
density = NULL, col.contour = "darkred", hex.bins = 50,
hex.pal = "viridis" ) +
coord_fixed(ratio = 1) # fixer le ratio/proportionRemarque : il est (très) utile d’avoir recours à coord_fixed(ratio = 1)afin de fixer le ratio des axes des X et Y et de conserver les proportions.
# nuage de TOUTES les modalités actives
# toujours dans le plan 1-2
# toutes les modalités : points = "all"
ggcloud_variables(mca_duval, axes = c(1,2), points = "all",
min.ctr = NULL, max.pval = 0.01, face = "pp",
shapes = TRUE, prop = NULL, textsize = 3, shapesize = 3,
col = NULL, col.by.group = TRUE, alpha = 1,
segment.alpha = 0.5, vlab = FALSE, legend = "right",
force = 1, max.overlaps = Inf) +
coord_fixed(ratio = 1)# Nuage des modalités actives qui contribuent à l'axe horizontal
# toujours dans le plan 1-2
# toutes les modalités : points = "besth"
ggcloud_variables(mca_duval, axes = c(1,2), points = "besth",
min.ctr = NULL, max.pval = 0.01, face = "pp",
shapes = TRUE, prop = NULL, textsize = 3, shapesize = 3,
col = NULL, col.by.group = TRUE, alpha = 1,
segment.alpha = 0.5, vlab = FALSE, legend = "right",
force = 1, max.overlaps = Inf)+
coord_fixed(ratio = 1)# Nuage des modalités actives qui contribuent à l'axe vertical
# toujours dans le plan 1-2
# toutes les modalités : points = "bestv"
ggcloud_variables(mca_duval, axes = c(1,2), points = "bestv",
min.ctr = NULL, max.pval = 0.01, face = "pp",
shapes = TRUE, prop = NULL, textsize = 3, shapesize = 3,
col = NULL, col.by.group = TRUE, alpha = 1,
segment.alpha = 0.5, vlab = FALSE, legend = "right",
force = 1, max.overlaps = Inf)+
coord_fixed(ratio = 1)# Nuage des modalités actives qui contribuent aux 2 axes
# toujours dans le plan 1-2
# toutes les modalités : points = "besthv"
ggcloud_variables(mca_duval, axes = c(1,2), points = "besthv",
min.ctr = NULL, max.pval = 0.01, face = "pp",
shapes = TRUE, prop = NULL, textsize = 3, shapesize = 3,
col = NULL, col.by.group = TRUE, alpha = 1,
segment.alpha = 0.5, vlab = FALSE, legend = "right",
force = 1, max.overlaps = Inf)+
coord_fixed(ratio = 1)# Superposition de variables supplémentaires 'au compte-gouttes'
# ici WW
# Step 1: on reprend le plot de base et on le stocke dans p
p <- ggcloud_variables(mca_duval, axes = c(1,2), points = "besthv",
min.ctr = NULL, max.pval = 0.01, face = "pp",
shapes = TRUE, prop = NULL, textsize = 3, shapesize = 3,
col = NULL, col.by.group = TRUE, alpha = 1,
segment.alpha = 0.5, vlab = FALSE, legend = "right",
force = 1, max.overlaps = Inf) +
coord_fixed(ratio = 1)
# Step 2: prépa ajout des variables supplémentaires
duval$WW <- as.factor(duval$WW)
freq(duval$WW) n % val%
Aucun 196 78.4 78.4
France 27 10.8 10.8
International 6 2.4 2.4
LesDeux 21 8.4 8.4
# ajout au plot existant
p <- ggadd_supvar(p, mca_duval, duval$WW, axes = c(1,2),
col = "black", shape = 1, prop = NULL,
textsize = 3, shapesize = 6,
segment = FALSE, vname = NULL)
# Step 3: imprimer le plot modifié
print(p)A VOUS DE JOUER !!